1 SETUP

1.1 Open libraries

library(gtsummary)
library(dplyr)
library(tidyverse)
library(survival)
library(survminer)
library(lubridate)
library(readxl)
library(janitor)
library(ggplot2)
library(tidycmprsk)
library(ggsurvfit)
library(tidymodels) 
library(sjPlot)
library(patchwork)
library(bstfun)
library(ComplexUpset)

library(writexl)
library(rms)

library(timeROC)
library(rms)
library(pec)
library(prodlim)
library(riskRegression)

library(plotly)
library(ggtext)

library(WeightIt)
library(cobalt)
library(survey)

2 DATA WRANGLING

2.1 Read in MM datasets

db_Kansas = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 2.11.2025 - due 5-1-25 (de-identified).xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  filter(!is.na(age_at_bs_ab_initiation)) %>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    alc_day_8 = as.double(alc_day_8),
    alc_day_15 = as.double(alc_day_15),
    alc_day_22 = as.double(alc_day_22),
    alc_c2d1 = as.double(alc_c2d1),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_Kansas, "Kansas_elra_data.xlsx")

db_CC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/CC bsAb Consortium Submission CCF 6-11-25.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    alc_day_8 = as.double(alc_day_8),
    alc_day_15 = as.double(alc_day_15),
    alc_day_22 = as.double(alc_day_22),
    alc_c2d1 = as.double(alc_c2d1),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_CC, "CC_elra_data.xlsx")

db_HUMC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/BsAb Consortium Data HUMC.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    alc_day_8 = as.double(alc_day_8),
    alc_day_15 = as.double(alc_day_15),
    alc_day_22 = as.double(alc_day_22),
    alc_c2d1 = as.double(alc_c2d1),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_HUMC, "HUMC_elra_data.xlsx")

db_MCC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data MCC-deidentified.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )


db_Huntsman = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/Utah_Huntsman_bsAb Consortium Data FINAL_5.4.25.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc)/1000,
    baseline_alc = as.double(baseline_alc)/1000,
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_Huntsman, "Huntsman_elra_data.xlsx")


db_FHCC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/FHCC Bispecifics (5.5.2025).xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    date_of_infection_80 = ymd(date_of_infection_80),
    baseline_anc = as.double(baseline_anc),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )


db_UTSW = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 2.11.2025 1.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3")) %>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles)
  )

# Write to Excel
#write_xlsx(db_UTSW, "UTSW_elra_data.xlsx")


db_MDACC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/MDACC BITE ASH 2025.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1),
    total_number_paraskelatal_plasmacytomas = as.double(total_number_paraskelatal_plasmacytomas)
  )

# Write to Excel
#write_xlsx(db_MDACC, "MDACC_elra_data.xlsx")


db_Roswell = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 5.8.2025_NN.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1),
    total_number_paraskelatal_plasmacytomas = as.double(total_number_paraskelatal_plasmacytomas)
  )

# Write to Excel
#write_xlsx(db_Roswell, "Roswell_elra_data.xlsx")

# Combine all into a single dataframe
combined_db <- bind_rows(
  list(
    HUMC = db_HUMC,
    MCC = db_MCC,
    Huntsman = db_Huntsman,
    FHCC = db_FHCC,
    UTSW = db_UTSW,
    MDACC = db_MDACC,
    Roswell = db_Roswell,
    CC = db_CC,
    Kansas = db_Kansas
  ),
  .id = "site"
) %>%
  mutate(
    prior_bcma_directed_therapy_yes_no = ifelse(!is.na(most_recent_bcma_agent_name),1,0),
    bs_ab_name = case_when(
      bs_ab_name == 1 ~ "tec",
      bs_ab_name == 3 ~ "elra",
      .default = NULL
    ),
    max_crs_grade_0_5 = ifelse(is.na(max_crs_grade_0_5), 0,max_crs_grade_0_5), ## Set NAs to 0
    max_icans_grade = ifelse(is.na(max_icans_grade), 0, max_icans_grade ) ## Set NAs to 0
  )

# Convert all datetime columns to Date
combined_db_clean <- combined_db %>%
  mutate(across(where(lubridate::is.POSIXt), as.Date))

# Write to Excel
write_xlsx(combined_db_clean, "elra_tec_data.xlsx")


# Write to Excel
# write_xlsx(combined_db_clean %>% filter(prior_bcma_directed_therapy_yes_no != prior_bcma_directed_therapy), "elra_data_discordant.xlsx")

## Several institutions do not have any cases of elranatamab

# db_UABMC = read_excel(path = "Data/bsAb Consortium Data 5.8.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_Stanford = read_excel(path = "Data/De-identified Stanford bsAb Consortium Data Template 4.30.2025.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_LCI = read_excel(path = "Data/bsAb Consortium Data LCI.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_COH = read_excel(path = "Data/bsAb Consortium Data Template 2.11.2025_COHv2.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_MUSC = read_excel(path = "Data/USMM bsAb Consortium Data Template 2.11.2025.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_Mayo = read_excel(path = "Data/Bispecific_consortium_20250521 NoPHI.xlsx", sheet = "bsAb_DATA") %>%
#   filter(`bsAb name` == 3)%>%
#   clean_names()

2.2 Calculate CAR-HEMATOTOX score and risk level

combined_db <- combined_db %>% 
  mutate(
    # baseline_anc = as.numeric(baseline_anc),
    # baseline_hgb,
    # baseline_pl_ts,
    # baseline_crp_prior_to_bs_ab_initiation,
    # baseline_ferritin_prior_to_bs_ab_initiation,
    HT_Score = 
      (baseline_anc <= 1.2) +
      (baseline_hgb <= 9.0) +
      (baseline_pl_ts >= 76 & baseline_pl_ts <= 175) +
      (baseline_crp_prior_to_bs_ab_initiation >= 3.0) +
      (baseline_ferritin_prior_to_bs_ab_initiation >= 650 & baseline_ferritin_prior_to_bs_ab_initiation <= 2000) +
      (baseline_pl_ts <= 75) + ## Add an extra point for Plt <=75
      (baseline_ferritin_prior_to_bs_ab_initiation > 2000), ## Add an extra point for ferritin >2000
    
    # Classify the risk based on the score
    HT_Risk = factor(
      case_when(
        HT_Score >= 2 ~ "HThigh",
        baseline_pl_ts <=75 ~ "HThigh", ## Equivalent to at least 2 points
        baseline_ferritin_prior_to_bs_ab_initiation > 2000 ~ "HThigh", ## Equivalent to at least 2 points
        (baseline_anc <= 1.2) + (baseline_hgb <= 9.0) + (baseline_pl_ts >= 76 & baseline_pl_ts <= 175) >=2 ~ "HThigh",
        HT_Score <2 ~ "HTlow",
        .default = "NA"
      ),
      levels = c("HTlow", "HThigh")
    )
)
write_xlsx(combined_db %>% select(site, patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate, baseline_anc, baseline_hgb, baseline_pl_ts, baseline_crp_prior_to_bs_ab_initiation, baseline_ferritin_prior_to_bs_ab_initiation, HT_Score, HT_Risk), "HT_score.xlsx")

2.3 Identify severe infections

combined_db <- combined_db %>% 
  mutate(
    severe_infection = ifelse(infection_severity_73 %in% 1 | infection_severity_76 %in% 1 | infection_severity_79 %in% 1| infection_severity_82 %in% 1,1,0),
    first_severe_infection_date = case_when(
      infection_severity_73 %in% 1 ~ date_of_infection_74,
      infection_severity_76 %in% 1 ~ date_of_infection_77,
      infection_severity_79 %in% 1 ~ date_of_infection_80,
      infection_severity_82 %in% 1 ~ date_of_infection_83,
      .default = NULL
    ),
    first_severe_infection_type = case_when(
      infection_severity_73 %in% 1 ~ number_infection_1_bacterial_viral_fungal,
      infection_severity_76 %in% 1 ~ number_infection_2_bacterial_viral_fungal,
      infection_severity_79 %in% 1 ~ number_infection_3_bacterial_viral_fungal,
      infection_severity_82 %in% 1 ~ number_infection_4_bacterial_viral_fungal,
      .default = NULL
    )
    
    
  )

3 PROPENSITY SCORES

3.1 Regression analysis for treatment allocation

theme_gtsummary_compact()

## Setting theme "Compact"

variables <- c("elra", "Age", "Female sex", "ECOG 2+", "EMD", "Hgb", "LDH", "Prior BCMA")

data1 <- combined_db %>% 
  drop_na(age_at_bs_ab_initiation, ecog_ps_at_bs_ab_initiation_0_5, just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk, baseline_hgb, baseline_ldh_prior_to_bs_ab_initiation) %>%
  transmute(
    Age = age_at_bs_ab_initiation,
    "Female sex" = sex_0_male_1_female,
    elra = ifelse(bs_ab_name == "elra", 1, 0),
    "ECOG 2+" = ifelse(ecog_ps_at_bs_ab_initiation_0_5 >=2,1,0),
    EMD = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,1,0),
    Hgb = as.double(baseline_hgb),
    LDH = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Prior BCMA" = prior_bcma_directed_therapy_yes_no
  )

uv_tab <- tbl_uvregression(
    data1[c(variables)],
    method = glm,
    y = "elra",
    method.args = list(family = binomial),
    exponentiate = TRUE
  ) %>%
  sort_p()


mv_tab<-glm(elra ~ Age + `Female sex` + `ECOG 2+` + EMD + Hgb + LDH + `Prior BCMA`, data=data1, family=binomial) %>%
  tbl_regression(exponentiate=TRUE)

tbl_merge(list(uv_tab, mv_tab), tab_spanner = c("**Univariate**", "**Multivariable**"))

Characteristic	Univariate				Multivariable
Characteristic	N	OR¹	95% CI¹	p-value	OR¹	95% CI¹	p-value
Hgb	371	1.12	1.01, 1.25	0.029	1.14	1.02, 1.28	0.026
EMD	371	1.69	0.95, 2.96	0.069	1.71	0.96, 3.04	0.068
Age	371	1.02	1.00, 1.04	0.087	1.01	0.99, 1.04	0.3
Female sex	371	1.43	0.93, 2.21	0.11	1.43	0.92, 2.24	0.11
Prior BCMA	371	0.76	0.49, 1.17	0.2	0.78	0.50, 1.23	0.3
ECOG 2+	371	1.19	0.75, 1.88	0.4	1.27	0.78, 2.07	0.3
LDH	371	1.00	1.00, 1.00	0.8	1.00	1.00, 1.00	>0.9
¹ OR = Odds Ratio, CI = Confidence Interval

3.2 IPTW

IPTW.data <- combined_db %>% 
  drop_na(age_at_bs_ab_initiation, ecog_ps_at_bs_ab_initiation_0_5, just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk, baseline_hgb, baseline_ldh_prior_to_bs_ab_initiation) %>%
  mutate(
    Age = age_at_bs_ab_initiation,
    "Female sex" = sex_0_male_1_female,
    bsAb = bs_ab_name,
    "Received full dose" = ifelse(!is.na(date_of_first_bs_ab_full_dose),1,0),
    "ECOG" = case_when(
      ecog_ps_at_bs_ab_initiation_0_5 %in% c("0","1") ~"0-1",
      ecog_ps_at_bs_ab_initiation_0_5 ==2 ~"2",
      ecog_ps_at_bs_ab_initiation_0_5 >=3 ~"3+",
      .default = NULL),
    EMD = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,1,0),
    Hgb = as.double(baseline_hgb),
    LDH = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Prior BCMA" = prior_bcma_directed_therapy_yes_no
  )

## Use the WeightIt function to generate propensity scores and weights
W.out <- weightit(bsAb ~ Age + `ECOG` + EMD + Hgb + LDH + `Prior BCMA`,
                  data=IPTW.data, 
                  method = "glm",    ## Methods: glm, gbm, cbps, npcbps, ebal, optweight, super, bart, energy
                  estimand = "ATE",
                  stabilize = TRUE)

bal.tab(W.out, stats = c("m", "v"), thresholds = c(m=0.25), binary = "std")

## Balance Measures
##                Type Diff.Adj     M.Threshold V.Ratio.Adj
## prop.score Distance   0.0050 Balanced, <0.25      0.9733
## Age         Contin.   0.0090 Balanced, <0.25      0.9595
## ECOG_0-1     Binary  -0.0121 Balanced, <0.25           .
## ECOG_2       Binary   0.0152 Balanced, <0.25           .
## ECOG_3+      Binary  -0.0027 Balanced, <0.25           .
## EMD          Binary  -0.0034 Balanced, <0.25           .
## Hgb         Contin.  -0.0186 Balanced, <0.25      1.0974
## LDH         Contin.  -0.0276 Balanced, <0.25      1.0961
## Prior BCMA   Binary  -0.0119 Balanced, <0.25           .
## 
## Balance tally for mean differences
##                     count
## Balanced, <0.25         9
## Not Balanced, >0.25     0
## 
## Variable with the greatest mean difference
##  Variable Diff.Adj     M.Threshold
##       LDH  -0.0276 Balanced, <0.25
## 
## Effective sample sizes
##              elra    tec
## Unadjusted 123.   248.  
## Adjusted   113.38 242.39

summary(W.out)

##                   Summary of weights
## 
## - Weight ranges:
## 
##         Min                                  Max
## elra 0.4745 |---------------------------| 2.3134
## tec  0.7811      |-------------|          1.7418
## 
## - Units with the 5 most extreme weights by group:
##                                         
##          123    102     81     91      9
##  elra 1.5458 1.5896  1.653 1.7607 2.3134
##          169     71    274    258    271
##   tec 1.4257 1.4453 1.5215 1.6783 1.7418
## 
## - Weight statistics:
## 
##      Coef of Var   MAD Entropy # Zeros
## elra       0.292 0.222   0.040       0
## tec        0.152 0.115   0.011       0
## 
## - Effective Sample Sizes:
## 
##              elra    tec
## Unweighted 123.   248.  
## Weighted   113.38 242.39

love_plot <- love.plot(W.out,
          binary = "std",
          line = TRUE, 
          abs = TRUE,
          thresholds = c(m=0.2),
          sample.names = c("Unweighted", "IPTW"),
          title = "A)",
          var.order = "unadjusted",
            var.names = c(
              prop.score = "Distance"
          ))+ theme(plot.title = element_text(hjust = 0))
love_plot

weight_df = data.frame(W.out$weights, W.out$treat)

IPTW.data$weights <- W.out$weights
  
IPTW.data$ps <- W.out$ps


write.csv(data1, "IPTW.dataset.csv")

3.3 Scatterplot of weights

Weights_plot <- ggplot(
  IPTW.data, aes(x=ps, y = weights, color = bsAb)) + 
  geom_point(shape = 16) +
  theme_classic() + 
  scale_y_continuous(breaks = seq(0,4 , by = 0.5)) +
  labs(x = "Propensity score", y = "Weights") + 
  ggtitle("B)")
Weights_plot

4 TABLES

4.1 Site contributions

theme_gtsummary_compact()

combined_db %>%
  transmute(
    "Site" = site,
    bs_ab_name,
  ) %>%
  tbl_summary(
    by = bs_ab_name,
    missing = "no",
    sort = list(all_categorical() ~ "frequency"),
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)"
    )
  ) %>%
  bold_labels()

Characteristic	elra N = 130¹	tec N = 277¹
Site
MCC	53 (41%)	72 (26%)
MDACC	12 (9.2%)	57 (21%)
CC	20 (15%)	47 (17%)
UTSW	10 (7.7%)	40 (14%)
FHCC	7 (5.4%)	24 (8.7%)
Huntsman	9 (6.9%)	22 (7.9%)
Kansas	4 (3.1%)	13 (4.7%)
HUMC	10 (7.7%)	1 (0.4%)
Roswell	5 (3.8%)	1 (0.4%)
¹ n (%)

4.2 Baseline characteristics

theme_gtsummary_compact()

combined_db %>%
  transmute(
    #"Site" = site,
    "Age (years)" = age_at_bs_ab_initiation,
    "Female sex" = sex_0_male_1_female == 1,
    "Race" = factor(case_when(
      is.na(race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5) ~ NA,
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 0 ~ "White",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 1 ~ "Black",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 2 ~ "AAPI",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 3 ~ "American Indian",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 4 ~ "Other",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 5 ~ NA
    ), levels = c("White","Black","AAPI","American Indian","Other",NA)),
    #ecog_ps_at_bs_ab_initiation_0_5,
    "ECOG" = case_when(
      ecog_ps_at_bs_ab_initiation_0_5 == 0 ~ "0",
      ecog_ps_at_bs_ab_initiation_0_5 == 1 ~ "1",
      ecog_ps_at_bs_ab_initiation_0_5 == 2 ~ "2",
      ecog_ps_at_bs_ab_initiation_0_5 >= 3 ~ "3+"
    ),
    "LDH (U/L)" = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Hgb (g/dL)" = baseline_hgb,
    "Heavy chain" = factor(case_when(
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgG" ~ "IgG",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgM" ~ "IgM",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgA" ~ "IgA",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgD" ~ "IgD",
      .default = "None"
      ), levels = c("IgG", "IgA","IgM","IgD","None")),
    "True EMD" = just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,
    "Oligo or non-secretory" = oligo_or_non_secretory_at_bs_ab_initiation_0_no_1_yes,
    "≥1 HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes==1 | t_4_14_prior_to_bs_ab_0_no_1_yes == 1 | t_14_16_prior_to_bs_ab_0_no_1_yes == 1 | amp_1_q_only_0_no_1_yes == 1 | gain_amp1q21_prior_to_bs_ab_0_no_1_yes == 1,1,0),
    "del(17p)" = deletion_17p_prior_to_bs_ab_0_no_1_yes,
    "t(4;14)" = t_4_14_prior_to_bs_ab_0_no_1_yes,
    "t(14;16)" = t_14_16_prior_to_bs_ab_0_no_1_yes,
    "gain/amp(1q)" = gain_amp1q21_prior_to_bs_ab_0_no_1_yes,
    "amp(1q)" = amp_1_q_only_0_no_1_yes,
    "Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
    "Prior ASCT" = prior_auto_sct_no_0_yes_1,
    "Prior CAR-T" = prior_cart_0_no_1_yes,
    "Bortezomib refractory" = bortezomib_status == 1,
    "Lenalidomide refractory" = lenalidomide_status == 1,
    "Pomalidomide refractory" = pomalidomide_status == 1,
    "Anti-CD38 refractory" = daratumumab_or_isatuximab_status == 1,
    "Triple refractory" = triple_class_refractory_i_mi_d_pi_anti_cd38 == 1,
    "Penta refractory" = penta_refractory_2_p_is_2_i_mi_ds_anti_cd38 == 1,
    "Prior BCMA-directed therapy" = !is.na(most_recent_bcma_agent_name),
    "Most recent BCMA" = factor(case_when(
      most_recent_bcma_agent_name %in% c("Abecma", "ABECMA","Abecma on trial","Ide-cel") ~ "CAR-T",
      most_recent_bcma_agent_name %in% c("Carvykti","carvykti", "Carvikty", "Clita-cel", "Cilta-cel") ~ "CAR-T",
      most_recent_bcma_agent_name %in% c("CC98633 CAR-T","ddBCMA") ~ "CAR-T",
      most_recent_bcma_agent_name %in% c("BLENREP", "blenrep", "Belantamab") ~ "ADC",
      most_recent_bcma_agent_name %in% c("Teclistamab") ~ "TCE"
    ), levels = c("CAR-T","TCE", "ADC")),
    "Trial eligible" = eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,
    "Received full dose" = ifelse(!is.na(date_of_first_bs_ab_full_dose),1,0),
     "Baseline ferritin (ng/mL)" = as.double(baseline_ferritin_prior_to_bs_ab_initiation),
     "Baseline Hgb" = as.double(baseline_hgb),
     "Baseline LDH" = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Duration of prior BCMA" = as.numeric(ymd(date_of_last_exposure_to_prior_bcma_agent) - ymd(date_of_first_exposure_to_prior_bcma_agent)),
    bs_ab_name
  ) %>%
  tbl_summary(
    by = bs_ab_name,
    missing = "no",
    #sort = list(all_categorical() ~ "frequency"),
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)",
      "Age (years)" ~ "{median} ({min} to {max})"
    )
  ) %>%
  add_p() %>%
  bold_labels() %>%
  add_variable_grouping(
    "High-risk cytogenetic abnormality at any time" = c("del(17p)","t(4;14)", "t(14;16)", "gain/amp(1q)", "amp(1q)"),
    "Refractory status" = c("Bortezomib refractory", "Lenalidomide refractory", "Pomalidomide refractory", "Anti-CD38 refractory")
  )

Characteristic	elra N = 130¹	tec N = 277¹	p-value²
Age (years)	71 (39 to 95)	69 (35 to 88)	0.085
Female sex	74 (57%)	128 (46%)	0.044
Race			0.14
White	102 (82%)	211 (76%)
Black	16 (13%)	54 (20%)
AAPI	3 (2.4%)	7 (2.5%)
American Indian	2 (1.6%)	0 (0%)
Other	2 (1.6%)	4 (1.4%)
ECOG			0.5
0	22 (17%)	52 (20%)
1	60 (48%)	130 (49%)
2	30 (24%)	65 (25%)
3+	14 (11%)	18 (6.8%)
LDH (U/L)	223 (177, 277)	208 (166, 287)	0.4
Hgb (g/dL)	10.05 (8.70, 11.70)	9.60 (8.30, 11.30)	0.020
Heavy chain			<0.001
IgG	75 (58%)	122 (44%)
IgA	28 (22%)	47 (17%)
IgM	2 (1.5%)	1 (0.4%)
IgD	2 (1.5%)	1 (0.4%)
None	23 (18%)	106 (38%)
True EMD	28 (22%)	36 (13%)	0.028
Oligo or non-secretory			0.3
0	116 (91%)	259 (94%)
1	12 (9.4%)	15 (5.4%)
2	0 (0%)	2 (0.7%)
≥1 HRCA	85 (69%)	164 (63%)	0.2
High-risk cytogenetic abnormality at any time
del(17p)	31 (25%)	70 (27%)	0.6
t(4;14)	17 (13%)	34 (13%)	>0.9
t(14;16)	8 (6.3%)	15 (5.8%)	0.8
gain/amp(1q)	70 (56%)	128 (49%)	0.2
amp(1q)	9 (7.9%)	18 (6.9%)	0.7
Prior LOTs	6.00 (4.00, 8.00)	6.00 (4.00, 8.00)	0.8
Prior ASCT	66 (51%)	163 (59%)	0.13
Prior CAR-T	54 (42%)	103 (37%)	0.4
Refractory status
Bortezomib refractory	87 (67%)	226 (82%)	<0.001
Lenalidomide refractory	110 (85%)	238 (86%)	0.7
Pomalidomide refractory	96 (74%)	226 (82%)	0.062
Anti-CD38 refractory	116 (89%)	267 (97%)	0.002
Triple refractory	118 (91%)	253 (92%)	0.8
Penta refractory	64 (49%)	141 (51%)	0.7
Prior BCMA-directed therapy	64 (49%)	152 (55%)	0.3
Most recent BCMA			0.012
CAR-T	54 (84%)	57 (70%)
TCE	6 (9.4%)	5 (6.2%)
ADC	4 (6.3%)	19 (23%)
Trial eligible	28 (22%)	63 (23%)	0.8
Received full dose	121 (93%)	268 (97%)	0.093
Baseline ferritin (ng/mL)	293 (132, 832)	393 (131, 1,254)	0.3
Baseline Hgb	10.05 (8.70, 11.70)	9.60 (8.30, 11.30)	0.020
Baseline LDH	223 (177, 277)	208 (166, 287)	0.4
Duration of prior BCMA	0 (0, 22)	0 (0, 126)	0.034
¹ Median (Min to Max); n (%); Median (Q1, Q3)
² Wilcoxon rank sum test; Pearson’s Chi-squared test; Fisher’s exact test

4.3 Unweighted and IPTW baseline characteristics

theme_gtsummary_compact()

## Setting theme "Compact"

data1 = IPTW.data %>%
  transmute(
    "Age (years)" = age_at_bs_ab_initiation,
    "Female sex" = sex_0_male_1_female == 1,
    "Race" = factor(case_when(
      is.na(race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5) ~ NA,
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 0 ~ "White",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 1 ~ "Black",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 2 ~ "AAPI",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 3 ~ "American Indian",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 4 ~ "Other",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 5 ~ NA
    ), levels = c("White","Black","AAPI","American Indian","Other",NA)),
    #ecog_ps_at_bs_ab_initiation_0_5,
    "ECOG" = case_when(
      ecog_ps_at_bs_ab_initiation_0_5 == 0 ~ "0",
      ecog_ps_at_bs_ab_initiation_0_5 == 1 ~ "1",
      ecog_ps_at_bs_ab_initiation_0_5 == 2 ~ "2",
      ecog_ps_at_bs_ab_initiation_0_5 >= 3 ~ "3+"
    ),
    "LDH (U/L)" = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Hgb (g/dL)" = baseline_hgb,
    "Baseline ferritin (ng/mL)" = as.double(baseline_ferritin_prior_to_bs_ab_initiation),
    "Heavy chain" = factor(case_when(
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgG" ~ "IgG",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgM" ~ "IgM",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgA" ~ "IgA",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgD" ~ "IgD",
      .default = "None"
      ), levels = c("IgG", "IgA","IgM","IgD","None")),
    "True EMD" = just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,
    "≥1 HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes==1 | t_4_14_prior_to_bs_ab_0_no_1_yes == 1 | t_14_16_prior_to_bs_ab_0_no_1_yes == 1 | amp_1_q_only_0_no_1_yes == 1 | gain_amp1q21_prior_to_bs_ab_0_no_1_yes == 1,1,0),
    "Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
    "Prior ASCT" = prior_auto_sct_no_0_yes_1,
    "Prior CAR-T" = prior_cart_0_no_1_yes,
    "Triple refractory" = triple_class_refractory_i_mi_d_pi_anti_cd38 == 1,
    "Penta refractory" = penta_refractory_2_p_is_2_i_mi_ds_anti_cd38 == 1,
    "Prior BCMA-directed therapy" = !is.na(most_recent_bcma_agent_name),
    "Trial eligible" = eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,
    "Received full dose" = ifelse(!is.na(date_of_first_bs_ab_full_dose),1,0),
    bs_ab_name,
    weights
  )

table_unweighted <- data1 %>% select(-weights) %>%
  tbl_summary(
    by = bs_ab_name,
    missing = "no",
    #sort = list(all_categorical() ~ "frequency"),
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)",
      "Age (years)" ~ "{median} ({min} to {max})"
    )
  ) %>%
  add_p() %>%
  bold_labels()

table_IPTW <- survey::svydesign(
    id = ~0,
    weights = ~weights,
    data = data1,
    fpc = NULL
  ) %>%
  tbl_svysummary(
    by = bs_ab_name,
    missing = "no",
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)"),
    include = c("Age (years)", "Female sex","Race","ECOG","LDH (U/L)","Hgb (g/dL)","Baseline ferritin (ng/mL)", "Heavy chain","True EMD","≥1 HRCA","Prior LOTs","Prior ASCT","Prior CAR-T","Triple refractory","Penta refractory","Trial eligible","Prior BCMA-directed therapy", "Received full dose")
  ) %>%
  add_p() %>%
    bold_labels()

## The following warnings were returned during `add_p()`:
## ! For variable `Heavy chain` (`bs_ab_name`) and "statistic" and "p.value"
##   statistics: Chi-squared approximation may be incorrect
## ! For variable `Race` (`bs_ab_name`) and "statistic" and "p.value" statistics:
##   Chi-squared approximation may be incorrect

tbl_merge(list(table_IPTW, table_unweighted), tab_spanner = c("**IPTW**", "**Unweighted**"))

Characteristic	IPTW			Unweighted
Characteristic	elra N = 123¹	tec N = 248¹	p-value²	elra N = 123³	tec N = 248³	p-value⁴
Age (years)	70 (64, 76)	69 (62, 76)	0.9	71 (39 to 95)	68 (35 to 88)	0.060
Female sex	66 (54%)	114 (46%)	0.2	68 (55%)	115 (46%)	0.11
Race			0.2			0.15
White	94 (79%)	187 (76%)		95 (81%)	187 (76%)
Black	17 (14%)	50 (20%)		16 (14%)	51 (21%)
AAPI	3 (2.5%)	6 (2.3%)		3 (2.5%)	5 (2.0%)
American Indian	2 (1.8%)	0 (0%)		2 (1.7%)	0 (0%)
Other	2 (2.0%)	4 (1.8%)		2 (1.7%)	4 (1.6%)
ECOG			>0.9			0.5
0	23 (18%)	49 (20%)		21 (17%)	50 (20%)
1	62 (50%)	119 (48%)		59 (48%)	121 (49%)
2	28 (23%)	59 (24%)		29 (24%)	60 (24%)
3+	10 (8.2%)	20 (8.2%)		14 (11%)	17 (6.9%)
LDH (U/L)	223 (178, 280)	208 (167, 290)	0.4	222 (175, 277)	208 (167, 291)	0.5
Hgb (g/dL)	9.90 (8.50, 11.40)	9.80 (8.40, 11.50)	0.9	10.00 (8.70, 11.70)	9.60 (8.25, 11.20)	0.027
Baseline ferritin (ng/mL)	346 (140, 1,250)	379 (125, 1,219)	>0.9	293 (137, 832)	393 (128, 1,254)	0.3
Heavy chain			0.001			<0.001
IgG	72 (58%)	109 (44%)		72 (59%)	107 (43%)
IgA	26 (21%)	43 (17%)		26 (21%)	42 (17%)
IgM	1 (0.7%)	1 (0.4%)		1 (0.8%)	1 (0.4%)
IgD	2 (1.5%)	0 (0%)		2 (1.6%)	0 (0%)
None	23 (19%)	95 (38%)		22 (18%)	98 (40%)
True EMD	20 (16%)	39 (16%)	>0.9	26 (21%)	34 (14%)	0.067
≥1 HRCA	81 (70%)	151 (65%)	0.3	81 (70%)	152 (65%)	0.4
Prior LOTs	6.00 (5.00, 8.00)	6.00 (4.00, 8.00)	0.4	6.00 (4.00, 8.00)	6.00 (4.00, 8.00)	>0.9
Prior ASCT	66 (54%)	141 (57%)	0.6	61 (50%)	144 (58%)	0.12
Prior CAR-T	53 (43%)	90 (36%)	0.2	50 (41%)	94 (38%)	0.6
Triple refractory	112 (91%)	224 (91%)	0.9	112 (91%)	226 (91%)	>0.9
Penta refractory	61 (49%)	123 (49%)	>0.9	61 (50%)	126 (51%)	0.8
Trial eligible	24 (19%)	55 (22%)	0.5	26 (21%)	53 (21%)	>0.9
Prior BCMA-directed therapy	67 (54%)	133 (54%)	>0.9	60 (49%)	138 (56%)	0.2
Received full dose	114 (93%)	241 (97%)	0.059	114 (93%)	241 (97%)	0.045
¹ Median (Q1, Q3); n (%)
² Design-based KruskalWallis test; Pearson’s X^2: Rao & Scott adjustment
³ Median (Min to Max); n (%); Median (Q1, Q3)
⁴ Wilcoxon rank sum test; Pearson’s Chi-squared test; Fisher’s exact test

5 STACKED BAR GRAPH

5.1 FIGURE 1A: Toxicity, IPTW

data1 = IPTW.data %>%
  transmute(
    BsAb = as.factor(bs_ab_name),
    CRS.any = ifelse( max_crs_grade_0_5 != 0, 1, 0),
    ICANS.any = ifelse( max_icans_grade != 0,1,0),
    
    CRS.two = ifelse( max_crs_grade_0_5 == 2 | max_crs_grade_0_5 == 3 | max_crs_grade_0_5 == 4, 1, 0),
    ICANS.two = ifelse( max_icans_grade ==2 | max_icans_grade == 3 | max_icans_grade == 4,1,0),
    
    CRS.three = ifelse( max_crs_grade_0_5 == 3 | max_crs_grade_0_5 == 4, 1, 0),
    ICANS.three = ifelse( max_icans_grade == 3 | max_icans_grade == 4,1,0),
    weights
  )


weighted_table <- survey::svydesign(
    id = ~0,
    weights = ~weights,
    data = data1,
    fpc = NULL
  )


combined_data <- rbind(
  IPTW.data %>%
    count(bs_ab_name, max_crs_grade_0_5, wt=weights) %>%
    group_by(bs_ab_name) %>%
    mutate(
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      total = sum(n[max_crs_grade_0_5 != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_crs_grade_0_5 != 0]),
      Grade = factor(max_crs_grade_0_5, levels = c(0, 4, 3, 2, 1)),
      Category = "CRS"
    ),
  IPTW.data %>%
    count(bs_ab_name, max_icans_grade, wt=weights) %>%
    group_by(bs_ab_name) %>%
    mutate(
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      total = sum(n[max_icans_grade != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_icans_grade != 0]),
      Grade = factor(max_icans_grade, levels = c(0, 4, 3, 2, 1)),
      Category = "ICANS"
    )
) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )


# Any grade toxicity
CRS.pval <- as.numeric(svychisq(~CRS.any+BsAb, weighted_table)$p.value)
ICANS.pval <- as.numeric(svychisq(~ICANS.any+BsAb, weighted_table)$p.value)

f_labels <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.pval, ICANS.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )

# G2+ toxicity
CRS.two.pval <- as.numeric(svychisq(~CRS.two+BsAb, weighted_table)$p.value)
ICANS.two.pval <- as.numeric(svychisq(~ICANS.two+BsAb, weighted_table)$p.value)

f_labels_two <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.two.pval, ICANS.two.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )

# G3+ toxicity
CRS.three.pval <- as.numeric(svychisq(~CRS.three+BsAb, weighted_table)$p.value)
ICANS.three.pval <- as.numeric(svychisq(~ICANS.three+BsAb, weighted_table)$p.value)

f_labels_three <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.three.pval, ICANS.three.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )




toxplot_IPTW <- ggplot(combined_data) +
  aes(x = bs_ab_name, y = pct, fill = Grade) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(Grade != 0)) +
    geom_text(aes(label = paste0(signif(pct,2), "% (", round(n,0) ,")")), 
            data = . %>% filter(Grade != 0 & Grade != 4),
            position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (", signif(total,2),")"), x = bs_ab_name, y = total_pct + 6, vjust = 0), 
             data = . %>% filter(Grade == 1),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_fill_brewer(palette = "Pastel1") +
  #guides(fill = guide_legend(reverse = TRUE)) +
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    #legend.position = "bottom"
    legend.position = "right"
  ) + 
  coord_cartesian(ylim = c(0, 115)) +
  geom_text( aes(label = paste0("Grade 1+, p=", signif(label,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels ) + 
  geom_text( aes(label = paste0("Grade 2+, p=", signif(label,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_two ) + 
  geom_text( aes(label = paste0("Grade 3+, p=", signif(label,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_three ) + 
  ggtitle("A) IPTW")

toxplot_IPTW

5.2 FIGURE 1B: Toxicity, unweighted

combined_data <- rbind(
  combined_db %>%
    count(bs_ab_name, max_crs_grade_0_5) %>%
    group_by(bs_ab_name) %>%
    mutate(
      total = sum(n[max_crs_grade_0_5 != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_crs_grade_0_5 != 0]),
      Grade = factor(max_crs_grade_0_5, levels = c(0, 4, 3, 2, 1)),
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      Category = "CRS"
    ),
  combined_db %>%
    count(bs_ab_name, max_icans_grade) %>%
    group_by(bs_ab_name) %>%
    mutate(
      total = sum(n[max_icans_grade != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_icans_grade != 0]),
      Grade = factor(max_icans_grade, levels = c(0, 4, 3, 2, 1)),
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      Category = "ICANS"
    )
) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
    
  )




## Any-grade toxicity
CRS.d.test <- data.frame(
  results = ifelse(combined_db$max_crs_grade_0_5 !=0,1,0) ,
  test = combined_db$bs_ab_name
)
CRS.pval <- chisq.test(CRS.d.test$results,CRS.d.test$test, correct = FALSE)$p.value

ICANS.d.test <- data.frame(
  results = ifelse(combined_db$max_icans_grade !=0,1,0) ,
  test = combined_db$bs_ab_name
)
ICANS.pval <- chisq.test(ICANS.d.test$results,ICANS.d.test$test, correct = FALSE)$p.value


f_labels <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.pval, ICANS.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS") )
  )


## G2+ toxicity
CRS.two.d.test <- data.frame(
  results = ifelse(combined_db$max_crs_grade_0_5 == 2 | combined_db$max_crs_grade_0_5 == 3 | combined_db$max_crs_grade_0_5 == 4,1,0),
  test = combined_db$bs_ab_name
)
CRS.two.pval <- chisq.test(CRS.two.d.test$results,CRS.two.d.test$test, correct = FALSE)$p.value

ICANS.two.d.test <- data.frame(
  results = ifelse(combined_db$max_icans_grade == 2 | combined_db$max_icans_grade == 3 | combined_db$max_icans_grade == 4,1,0),
  test = combined_db$bs_ab_name
)
ICANS.two.pval <- chisq.test(ICANS.two.d.test$results,ICANS.two.d.test$test, correct = FALSE)$p.value

f_labels_two <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.two.pval, ICANS.two.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS") )
  )

## G3+ toxicity
CRS.three.d.test <- data.frame(
  results = ifelse( combined_db$max_crs_grade_0_5 == 3 | combined_db$max_crs_grade_0_5 == 4,1,0),
  test = combined_db$bs_ab_name
)
CRS.three.pval <- fisher.test(CRS.three.d.test$results,CRS.three.d.test$test)$p.value

ICANS.three.d.test <- data.frame(
  results = ifelse(combined_db$max_icans_grade == 3 | combined_db$max_icans_grade == 4,1,0),
  test = combined_db$bs_ab_name
)
ICANS.three.pval <- chisq.test(ICANS.three.d.test$results,ICANS.three.d.test$test, correct = FALSE)$p.value

f_labels_three <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.three.pval, ICANS.three.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS") )
  )



toxplot_unweighted <- ggplot(combined_data) +
  aes(x = bs_ab_name, y = pct, fill = Grade) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(Grade != 0)) +
    geom_text(aes(label = paste0(signif(pct,2), "% (",n,")")), 
            data = . %>% filter(Grade != 0 & Grade != 4),
            position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (",total,")"), x = bs_ab_name, y = total_pct + 6, vjust = 0), 
             data = . %>% filter(Grade == 1),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 100)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "none"
  ) + 
  coord_cartesian(ylim = c(0, 114)) +
  scale_fill_brewer(palette = "Pastel1", direction = 1) +
  geom_text( aes(label = paste0("Grade 1+, p=", signif(label,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels ) + 
  geom_text( aes(label = paste0("Grade 2+, p=", signif(label,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_two ) + 
  geom_text( aes(label = paste0("Grade 3+, p=", signif(label,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_three ) + 
  ggtitle("B) Unweighted")

toxplot_unweighted

5.3 ***FIGURE 1: IPTW and unweighted toicity

FIGURE1 <- toxplot_IPTW / toxplot_unweighted

ggsave("svg/FIGURE1.svg", FIGURE1, device = "svg",
       width = 7,
       height = 11,
       units = c("in"))

ggsave("tiff/FIGURE1.tiff", FIGURE1, device = "tiff",
       width = 7,
       height = 11,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

ggsave("jpeg/FIGURE1.jpeg", FIGURE1, device = "jpeg",
       width = 7,
       height = 11,
       dpi = 600,
       units = "in")

knitr::include_graphics("jpeg/FIGURE1.jpeg")

5.4 FIGURE 2A: Response, IPTW

data1 = IPTW.data %>%
  filter(!is.na(best_response)) %>%
  transmute(
    BsAb = as.factor(bs_ab_name),
    best_response = factor(best_response, levels =  c("sCR", "CR","VGPR","PR","MR","SD","PD")),
    ORR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR",1,0),
    VGPR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR",1,0),
    CR = ifelse(best_response == "sCR" | best_response == "CR",1,0),
    weights
  )

weighted_table <- survey::svydesign(
    id = ~0,
    weights = ~weights,
    data = data1,
    fpc = NULL
  )



data1 <-   IPTW.data %>%
  filter(!is.na(best_response)) %>%
  count(bs_ab_name, best_response, wt=weights) %>%
  group_by(bs_ab_name) %>%
  mutate(
    total = sum(n[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
    pct = prop.table(n) * 100,
    total_pct = sum(pct[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
    best_response = factor(best_response, levels =  c("sCR", "CR","VGPR","PR","MR", "SD","PD")),
    bs_ab_name = as.factor(bs_ab_name),
    Category = "Response rate"
  )



Response.pval <- as.numeric(svychisq(~ORR+BsAb, weighted_table)$p.value)
CR.pval <- as.numeric(svychisq(~CR+BsAb, weighted_table)$p.value)
VGPR.pval <- as.numeric(svychisq(~VGPR+BsAb, weighted_table)$p.value)

f_labels <- data.frame(Category = c("Response rate"), label = c(Response.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_CR <- data.frame(Category = c("Response rate"), label = c(CR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_VGPR <- data.frame(Category = c("Response rate"), label = c(VGPR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

responseplot_IPTW <- ggplot(data1) +
  aes(x = bs_ab_name, y = pct, fill = best_response) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR")) +
    geom_text(aes(label = paste0(signif(pct,2), "% (", round(n,0) ,")")), 
            data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR"| best_response == "PR"),
            position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (", signif(total,2),")"), x = bs_ab_name, y = total_pct + 7, vjust = 0), 
             data = . %>% filter(best_response == "PR"),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_fill_brewer(
    palette = "Pastel1",
    direction = -1,
    breaks = rev(levels(data1$best_response))
  ) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    #legend.position = "bottom"
    legend.position = "right"
  ) + 
  coord_cartesian(ylim = c(0, 114)) +
  geom_text( aes(label = paste0("ORR, p=", signif(Response.pval,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels  ) + 
  geom_text( aes(label = paste0("≥VGPR, p=", signif(VGPR.pval,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_VGPR  ) + 
  geom_text( aes(label = paste0("≥CR, p=", signif(CR.pval,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_CR  ) + 
  labs(fill = "Best response")+ 
  ggtitle("A) IPTW")

responseplot_IPTW

5.5 FIGURE 2B: Response, unweighted

data1 <- combined_db %>%
  filter(!is.na(best_response)) %>%
  transmute(
    best_response = factor(
      best_response, 
      levels =  c("sCR", "CR","VGPR", "PR","MR","SD","PD")),
    ORR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" |best_response == "PR",1,0),
    bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
    Category = "Response rate"
  )

combined_data <- combined_db %>%
  filter(!is.na(best_response)) %>%
  count(bs_ab_name, best_response) %>%
  group_by(bs_ab_name) %>%
  transmute(
    bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
    best_response = factor(
      best_response, 
      levels =  c("sCR", "CR","VGPR", "PR","MR","SD","PD")),
    total = sum(n[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
    pct = prop.table(n) * 100,
    n,
    total_pct = sum(pct[best_response == "sCR" | best_response == "CR" |best_response == "VGPR" | best_response == "PR"]),
    Category = "Response rate"
  )

## ORR
Response.d.test <- data.frame(
  results = ifelse(data1$ORR,1,0),
  test = data1$bs_ab_name
)
Response.pval <- chisq.test(Response.d.test$results,Response.d.test$test, correct = FALSE)$p.value

## VGPR
VGPR.d.test <- data.frame(
  results = ifelse(data1$best_response == "CR" | data1$best_response == "VGPR",1,0),
  test = data1$bs_ab_name
)
VGPR.pval <- chisq.test(VGPR.d.test$results,VGPR.d.test$test, correct = FALSE)$p.value

## CR
CR.d.test <- data.frame(
  results = ifelse(data1$best_response == "CR",1,0),
  test = data1$bs_ab_name
)
CR.pval <- chisq.test(CR.d.test$results,CR.d.test$test, correct = FALSE)$p.value


f_labels <- data.frame(Category = c("Response rate"), label = c(Response.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_CR <- data.frame(Category = c("Response rate"), label = c(CR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_VGPR <- data.frame(Category = c("Response rate"), label = c(VGPR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )


responseplot_unweighted <- ggplot(combined_data) +
  aes(x = bs_ab_name, y = pct, fill = best_response) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(best_response == "sCR" | best_response == "CR" |best_response == "VGPR" | best_response == "PR")) +
  geom_text(aes(label = paste0(signif(pct,2), "% (",n,")")), 
             data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"),
             position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (",total,")"), x = bs_ab_name, y = total_pct + 7, vjust = 0), 
             data = . %>% filter(best_response == "PR"),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    #legend.position = "bottom"
    legend.position = "none"
  ) + 
  coord_cartesian(ylim = c(0, 114)) +
  scale_fill_brewer(palette = "Pastel1", direction = -1) +
  geom_text( aes(label = paste0("ORR, p=", signif(Response.pval,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels) + 
  geom_text( aes(label = paste0("≥VGPR, p=", signif(VGPR.pval,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_VGPR  ) + 
  geom_text( aes(label = paste0("≥CR, p=", signif(CR.pval,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_CR  ) + 
  labs(fill = "Best response")+ 
  ggtitle("B) Unweighted")

responseplot_unweighted

5.6 ***FIGURE 2: IPTW and unweighted response

FIGURE2 <- responseplot_IPTW / responseplot_unweighted

ggsave("svg/FIGURE2.svg", FIGURE2, device = "svg",
       width = 7,
       height = 11,
       units = c("in"))

ggsave("tiff/FIGURE2.tiff", FIGURE2, device = "tiff",
       width = 7,
       height = 11,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

ggsave("jpeg/FIGURE2.jpeg", FIGURE2, device = "jpeg",
       width = 7,
       height = 11,
       dpi = 600,
       units = "in")

knitr::include_graphics("jpeg/FIGURE2.jpeg")

6 SURVIVAL ANALYSIS

6.1 Length of follow-up (Reverse KM, months)

data1 <- IPTW.data %>%
  transmute(
    bs_ab_name,
    Death = ifelse(death_yes_no,1,0),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Reverse_death = ifelse(death_yes_no == 1, 0,1),
    weights
  )

reverse_km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1)
reverse_km_OS

## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 
##     1, data = data1, weights = weights)
## 
##    15 observations deleted due to missingness 
##      records   n events median 0.95LCL 0.95UCL
## [1,]     356 356    215   13.4    12.5    14.9

reverse_km_OS_elra <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1 %>% filter(bs_ab_name == "elra"))
reverse_km_OS_elra

## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 
##     1, data = data1 %>% filter(bs_ab_name == "elra"), weights = weights)
## 
##    6 observations deleted due to missingness 
##      records   n events median 0.95LCL 0.95UCL
## [1,]     117 117   74.5   7.47    6.03       9

reverse_km_OS_tec <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1 %>% filter(bs_ab_name == "tec"))
reverse_km_OS_tec

## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 
##     1, data = data1 %>% filter(bs_ab_name == "tec"), weights = weights)
## 
##    9 observations deleted due to missingness 
##      records   n events median 0.95LCL 0.95UCL
## [1,]     239 239    140   17.1    14.9    20.1

6.2 Stratified by BsAb

data1 <- IPTW.data %>%
  transmute(
    site,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    bs_ab_name,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )

data_dor <- IPTW.data %>%
  filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
  transmute(
    weights,
    bs_ab_name,
    day_30_response,
    day_90_response,
    best_response,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_of_bs_ab_step_up_d1,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    day_30_response,
    day_90_response,
    best_response,
    time_to_best_response = case_when(
      best_response == day_30_response ~ 30,
      best_response == day_90_response ~ 90,
      .default = NA
    ),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
  ) %>%
  filter(!is.na(time_to_best_response)) %>%
  mutate(
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1)- time_to_best_response,
  ) %>%
  filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)

km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS

## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, 
##     data = data1, weights = weights)
## 
##    15 observations deleted due to missingness 
##                 records   n events median 0.95LCL 0.95UCL
## bs_ab_name=elra     117 117   42.8     NA     8.6      NA
## bs_ab_name=tec      239 239   98.6   21.5    16.4      NA

med_time <- surv_median(km_OS)

cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))

OS <- ggsurvplot(km_OS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Overall survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "A)",
           size = 0.5
) 

OS$plot <- OS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

OS

km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS

## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data1, weights = weights)
## 
##    12 observations deleted due to missingness 
##                 records   n events median 0.95LCL 0.95UCL
## bs_ab_name=elra     118 118   68.1   3.77    2.57    7.17
## bs_ab_name=tec      241 241  151.1   8.50    7.10   11.30

med_time <- surv_median(km_PFS)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))


PFS <- ggsurvplot(km_PFS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Progression-free survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "B)",
           size = 0.5
)

PFS$plot <- PFS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

PFS

km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR

## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data_dor, weights = weights)
## 
##    3 observations deleted due to missingness 
##                 records     n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      62  60.2   20.6   11.7    6.13      NA
## bs_ab_name=tec      118 120.7   56.1   13.7    9.27    22.3

med_time <- surv_median(km_DOR)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))


DOR <- ggsurvplot(km_DOR,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           #risk.table.fontsize = 4,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Duration of response",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "C)",
           size = 0.5
)

DOR$plot <- DOR$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

DOR

6.3 ***FIGURE 3: PFS/OS/DOR, stratified by trial elibility

FIGURE3 <- wrap_plots(
  (OS$plot / OS$table) + plot_layout(heights = c(4, 0.9)),
  (PFS$plot / PFS$table) + plot_layout(heights = c(4, 0.9)),
  (DOR$plot / DOR$table) + plot_layout(heights = c(4, 0.9)),
  ncol = 1  # Force vertical stacking
)

ggsave("svg/FIGURE3.svg", FIGURE3, device = "svg",
       width = 6.5,
       height = 8.5,
       units = c("in"))

ggsave("jpeg/FIGURE3.jpeg", FIGURE3, device = "jpeg",
       width = 6.5,
       height = 8.5,
       dpi = 600,
       units = c("in"))

ggsave("tiff/FIGURE3.tiff", FIGURE3, device = "tiff",
       width = 6.5,
       height = 8.5,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

knitr::include_graphics("jpeg/FIGURE3.jpeg")

6.4 Schoenfeld residuals and PH test

data1 <- IPTW.data %>%
  transmute(
    site,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    elra_ind = ifelse(bs_ab_name == "elra", 1, 0),  # 1 = elra, 0 = tec
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )

## ---- Cox models ----
cox_OS <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ elra_ind, weights = weights, data = data1)

cox_PFS <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ elra_ind, weights = weights, data = data1)

## ---- PH diagnostics ----
zph_OS  <- cox.zph(cox_OS)
zph_PFS <- cox.zph(cox_PFS)

# Global PH p-values (same as elra_ind right now)
p_OS  <- signif(zph_OS$table["GLOBAL",  "p"], 3)
p_PFS <- signif(zph_PFS$table["GLOBAL", "p"], 3)

par(mfrow = c(1, 2))

plot(
  zph_OS,
  main = paste0(
    "Schoenfeld residuals – OS\n",
    "Global PH test: p = ", p_OS
  )
)

plot(
  zph_PFS,
  main = paste0(
    "Schoenfeld residuals – PFS\n",
    "Global PH test: p = ", p_PFS
  )
)

par(mfrow = c(1, 1))

Elranatamab & Teclistamab IPTW Analysis

Andrew Portuguese

November 20, 2025

1 SETUP

1.1 Open libraries

2 DATA WRANGLING

2.1 Read in MM datasets

2.2 Calculate CAR-HEMATOTOX score and risk level

2.3 Identify severe infections

3 PROPENSITY SCORES

3.1 Regression analysis for treatment allocation

3.2 IPTW

3.3 Scatterplot of weights

4 TABLES

4.1 Site contributions

4.2 Baseline characteristics

4.3 Unweighted and IPTW baseline characteristics

5 STACKED BAR GRAPH

5.1 FIGURE 1A: Toxicity, IPTW

5.2 FIGURE 1B: Toxicity, unweighted

5.3 ***FIGURE 1: IPTW and unweighted toicity

5.4 FIGURE 2A: Response, IPTW

5.5 FIGURE 2B: Response, unweighted

5.6 ***FIGURE 2: IPTW and unweighted response

6 SURVIVAL ANALYSIS

6.1 Length of follow-up (Reverse KM, months)

6.2 Stratified by BsAb

6.3 ***FIGURE 3: PFS/OS/DOR, stratified by trial elibility

6.4 Schoenfeld residuals and PH test